%!TEX program = xelatex
% 完整编译: xelatex -> biber/bibtex -> xelatex -> xelatex
\documentclass[lang=cn,11pt,a4paper]{elegantpaper}
\usepackage{bm}

\title{LANT-Hybrid: 用于无碰撞等离子体模拟的混合粒子-流体代码}
\author{矫金龙 \\ 浙江大学，浙江近代物理中心，天文研究所}
%\institute{\href{https://elegantlatex.org/}{Elegant\LaTeX{} 项目组}}

\version{0.10}
\date{\zhtoday}


% 本文档命令
\usepackage{array}
\newcommand{\ccr}[1]{\makecell{{\color{#1}\rule{1cm}{1cm}}}}

\begin{document}

\maketitle

\begin{abstract}
本文为LANT-Hybrid代码的说明文档。
\keywords{混合模拟，无碰撞等离子体，离子-韦伯不稳定性}
\end{abstract}

\section{引言}



\section{控制方程（Basic Assumptions and Equations）}

这里采用的混合模拟方法是将离子作为粒子处理，电子作为保持准中性的流体来处理。首先需要求解离子的运动方程，即洛伦兹力方程（cgs单位制），
\begin{gather}
m_i\frac{d\bm{v}_p}{dt}=q_i(\bm{E}+\frac{1}{c}\bm{v}_p\times\bm{B}),\\
\frac{d\bm{r}_p}{dt}=\bm{v}_p.
\end{gather}
其中$\bm{v}_p$和$\bm{r}_p$分别是离子的速度和位置。电场强度$\bm{E}$和磁感应强度$\bm{B}$由ICS方程、法拉第定律和广义欧姆定律联合求出。这里磁场来自于两部分，即由离子-韦伯不稳定性产生的等离子体自生磁场和由变化电场感生出的感应磁场。这两部分磁场分别由ICS方程和法拉第定律描述。
\begin{gather}
\frac{c}{4\pi en_e}\nabla^2\bm{B}_w-\frac{e}{\eta m_e c}\bm{B}_w+\nabla\times\bm{u}=0,\\
(or)~\bm{B}_w=\frac{\eta m_e c}{e}\nabla\times\bm{u}, \\
\frac{\partial \bm{B}_f}{\partial t}=-c\nabla\times\bm{E}, \\
\bm{B}=\bm{B}_w+\bm{B}_f,\\
\bm{E}=-\frac{1}{en_e}\nabla p_e-\frac{1}{c}\bm{u}\times\bm{B}+\frac{1}{4\pi e n_e}(\nabla\times\bm{B})\times\bm{B}.
\end{gather}
其中，$\bm{B}_w$和$\bm{B}_f$分别是离子-韦伯不稳定性和法拉第电磁感应产生的磁场。
最后，为保持方程组封闭，还需要状态方程和等离子体准中性条件。状态方程采用绝热状态方程（其中$\Gamma$是绝热指数）。
\begin{gather}
p_e=p_{e0}(\frac{n_e}{n_{e0}})^{\Gamma}, \\
en_e=q_in_i.
\end{gather}
为了简化数值求解，需要将上述方程组进行无量纲化，下面是各个物理量的无量纲参数，
\begin{gather*}
\frac{\omega_{pe}}{c}\bm{r}_p=\bm{r}^*_p, \frac{1}{c}\bm{v}_p=\bm{v}^*_p, \omega_{pe}t=t^*, \frac{e}{m_e\omega_{pe}c}\bm{E}=\bm{E}^*, \frac{e}{m_e\omega_{pe}c}\bm{B}=\bm{B}^*, \\
\frac{4\pi e^2}{m_e\omega_{pe}^2}n_{e,i}=n^*_{e,i}, \frac{4\pi e^2}{m_e^2\omega_{pe}^2c^2}p_e=p_e^*, \frac{1}{m_e}m_i=m_i^*, \frac{1}{e}q_i=q_i^*.
\end{gather*}
利用上述无量纲参数，可以将上述方程组无量纲化为如下形式（为简便省略*号），
\begin{gather}
m_i\frac{d\bm{v}_p}{dt}=q_i(\bm{E}+\bm{v}_p\times\bm{B}),\\
\frac{d\bm{r}_p}{dt}=\bm{v}_p, \\
\frac{1}{n_e}\nabla^2\bm{B}_w-\frac{1}{\eta}\bm{B}_w+\nabla\times\bm{u}=0,\\
(or)~\bm{B}_w=\eta\nabla\times\bm{u}, \\
\frac{\partial\bm{B}_f}{\partial t}=-\nabla\times\bm{E}, \\
\bm{B}=\bm{B}_w+\bm{B}_f,\\
\bm{E}=-\frac{1}{n_e}\nabla p_e-\bm{u}\times\bm{B}+\frac{1}{n_e}(\nabla\times\bm{B})\times\bm{B}, \\
p_e=p_{e0}(\frac{n_e}{n_{e0}})^{\Gamma}, \\
n_e=q_in_i.
\end{gather}

\section{数值实现（Numerical Implementations）}

\subsection{差分格式}
粒子的位置和速度推动采用蛙跳格式（Leapfrog），即在时间步$N$，粒子位置$\bm{r}^N_p$、电场强度$\bm{E}^N$和磁感应强度$\bm{B}^N$已知，粒子速度$\bm{v}^{N-1/2}_p$在时间步$N-1/2$已知。在已知$\bm{E}^N$和$\bm{B}^N$的条件下，可以得到$\bm{r}^{N+1}_p$和$\bm{v}^{N+1/2}_p$。具体的粒子推动计算过程采用Boris提出的具有二阶精度（$O(\Delta t^2)$）的半加速-旋转-半加速算法。
\begin{gather}
\bm{v}^{-}=\bm{v}^{N-1/2}_p+\frac{q_i\Delta t}{2m_i}\bm{E}^N,\\
\bm{d}=\frac{q_i\Delta t}{2m_i}\bm{B}^N,\\
\bm{v}^{0}=\bm{v}^{-}+\bm{v}^{-}\times\bm{d},\\
\bm{s}=\frac{2\bm{d}}{1+\bm{d}^2},\\
\bm{v}^{+}=\bm{v}^{-}+\bm{v}^{0}\times\bm{s},\\
\bm{v}^{N+1/2}_p=\bm{v}^{+}+\frac{q_i\Delta t}{2m_i}\bm{E}^N,\\
\bm{r}^{N+1}_p=\bm{r}^N_p+\bm{v}^{N+1/2}_p\Delta t.
\end{gather}
对离子速度$\bm{v}^{N+1/2}_p$进行统计，可以得到等离子体离子流速$\bm{u}^{N+1/2}$。同样，对离子位置$\bm{r}^{N+1}_p$进行统计，可以得到等离子体离子密度$n^{N+1}_i$，进而可以得到电子密度$n^{N+1/2}_e=q_i(n^{N+1}_i+n^{N}_i)/2$。已知$\bm{u}^{N+1/2}$和$n^{N+1/2}_e$后，即可求解ICS方程，得到$\bm{B}_w^{N+1/2}$。ICS方程是二阶的椭圆型偏微分方程，可以采用迭代法高效求解。常用的迭代方法包括简单迭代、高斯-赛德尔（Gauss-Seidel）迭代以及逐次超松弛迭代等，考虑程序的并行化处理，这里采用简单迭代格式。
\begin{equation}
\begin{split}
\bm{B}_w^{m+1}(i,j,k)=\frac{\epsilon(i,j,k)}{\Delta x^2}[\bm{B}_w^{m}(i+1,j,k)+\bm{B}_w^{m}(i-1,j,k)]\\
+\frac{\epsilon(i,j,k)}{\Delta y^2}[\bm{B}_w^{m}(i,j+1,k)+\bm{B}_w^{m}(i,j-1,k)]\\
+\frac{\epsilon(i,j,k)}{\Delta z^2}[\bm{B}_w^{m}(i,j,k+1)+\bm{B}_w^{m}(i,j,k-1)]\\
+\epsilon(i,j,k)n_e(i,j,k)\nabla\times\bm{u}(i,j,k).
\end{split}
\end{equation}
其中系数$\frac{1}{\epsilon(i,j,k)}=\frac{2}{\Delta x^2}+\frac{2}{\Delta y^2}+\frac{2}{\Delta z^2}+\frac{n_e(i,j,k)}{\eta}$。迭代时，以$\bm{B}_w^{N-1/2}$的值作为迭代初值。

第$N+1$时间步的电场强度$\bm{E}^{N+1}$和磁感应强度$\bm{B}^{N+1}$的计算采用CAM-CL方法。具体如下，
\begin{gather}
\bm{B}_f^{N+1/2}=\bm{B}_f^N-\frac{\Delta t}{2}\nabla\times\bm{E}^N,\\
\bm{B}^{N+1/2}=\bm{B}_f^{N+1/2}+\bm{B}_w^{N+1/2},\\
\begin{split}
\bm{E}^{N+1/2}=-\frac{1}{n_e^{N+1/2}}\nabla p_e^{N+1/2}-\bm{u}^{N+1/2}\times\bm{B}^{N+1/2}\\
+\frac{1}{n_e^{N+1/2}}(\nabla\times\bm{B}^{N+1/2})\times\bm{B}^{N+1/2}\\
=\bm{F}(\bm{B}^{N+1/2},n_e^{N+1/2},\bm{u}^{N+1/2}),\\
\end{split}\\
\bm{B}_f^{N+1}=\bm{B}_f^{N+1/2}-\frac{\Delta t}{2}\nabla\times\bm{E}^{N+1/2},\\
\bm{B}^{N+1}=\bm{B}_f^{N+1}+\bm{B}_w^{N+1/2},\\
\bm{\tilde{E}}=\bm{F}(\bm{B}^{N+1},n_e^{N+1},\bm{u}^{N+1/2}),\\
\bm{\tilde{u}}=\bm{u}^{N+1/2}+\frac{q_i\Delta t}{2m_i}(\bm{\tilde{E}}+\bm{u}^{N+1/2}\times\bm{B}^{N+1}),\\
\bm{E}^{N+1}=\bm{F}(\bm{B}^{N+1},n_e^{N+1},\bm{\tilde{u}}).
\end{gather}

\subsection{2D3V差分格式}

\subsection{宏粒子}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{particle.jpg}
  \caption{宏粒子形状}
\end{figure}

\subsection{边界条件}

\subsection{并行策略}


\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{comm.jpg}
  \caption{并行区域划分与并行通信示意图}
\end{figure}



\section{结果与讨论}

\subsection{无碰撞静电激波}

\subsection{磁化无碰撞激波}

\subsection{等离子体剪切流磁场}

\subsection{离子-韦伯不稳定性磁场}

\section{总结}

\section{致谢}
This research was supported by the Fundamental Research Funds for the Central Universities (No. 226-2022-00216).

\nocite{*}
\printbibliography[heading=bibintoc, title=\ebibname]

\appendix
%\appendixpage
\addappheadtotoc

\end{document}
